A Simple and Efficient Regularization Method for 3D BEM: 
Application to Frequency-Domain Elastodynamics 



Patrick Dangla*, Jean-Frangois Semblat**0, Haihong Xiao**, Nicolas Dele-pine** 

* LCPC/LMSGC, 2 allee Kepler, 77420 Champs-sur-Marne, France 
" LCPC, 58 bd Lefebvre, 75732 Paris Cedex 15, France 

Abstract 

An efficient and easy-to-implement method is proposed to regularize integral equations in the 3D Bound- 
ary Element Method. The method takes advantage of an assumed three-noded triangle discretization of the 
boundary surfaces. The method is based on the derivation of analytical expressions of singular integrals. To 
demonstrate the accuracy of the method, three elastodynamic problems are numerically worked out in fre- 
quency domain: cavity under harmonic pressure, diffraction of a plane wave by a spherical cavity, amplification 
of seismic waves in a semi-spherical alluvial basin (the second one is also investigated in time domain). The 
numerical results are compared to (semi-) analytical solutions; a close agreement is found for all problems 
showing the very good accuracy of the proposed method. 
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Regularization of Boundary Integral Equations 

In contrast to other discretization methods, the Boun- 
dary Element Method involves singular integrals. Ac- 
cording to the power, there are three kinds of singu- 
larities that depend on whether integrability is defined 

(i) in the ordinary Riemann sense (weak singularity), 

(ii) in the Cauchy principal value sense (strong singu- 
larity) or (Hi) in the Hadamard finite part sense (hyper- 
singularity) . Strong singular integrals appear in the Or- 
dinary Boundary Integral Equations (OBIE) while the 
Derivative Boundary Integral Equations (DBIE) involve 
both the strong and hyper-singular integrals. Strong 
and hyper-singular integrals have to be converted to 
regular ones in the regularization of the BEM formu- 
lations (Tanaka et al, 1994; Sladek and Sladek, 1996, 
1998). Strictly speaking, weak singularity is not treated 
by regularization. However, from the point of view of 
numerical integrations, one should devote a great atten- 
tion to the evaluation of these integrals because stan- 
dard integration quadratures fail in accuracy (Lachat 
and Watson, 1976; Sladek et al, 1997, 2001; Manolis 
and Beskos, 1988). Therefore each type of singularity 
has to be treated by appropriate techniques. Most of the 
researches has dealt with strong (Bonnet and Bui, 1993) 
and hyper-singularities. Some methods have been pro- 
posed in the literature to treat these singular integrals 
(Sladek and Sladek, 1996; Niu and Zhou, 2004; Guig- 
giani and Gigante, 1990; Guiggiani et al, 1992; Chen 
and Hong, 1999; Bonnet, 1999; Bui et al, 1985; Bonnet 
and Xiao, 1995; Aubry and Clouteau, 1991; Xiao, 1994; 
Guiggiani, 1994). It is noticed that the regularization 
can be performed either before or after the discretiza- 
tion, i.e. in the global or local (intrinsic) coordinate 
space, as observed in some papers mentioned above. A 
comprehensive review of BEM in dynamic analysis has 
been proposed by Beskos (Beskos, 1997). 

In this paper, the regularization is performed in the 
global coordinate space after the discretization of the 
geometry. Herein, only strong and weak singularities 
of the ordinary boundary integral equations are dealt 
with. The method takes advantage of an assumed three- 
noded triangle element for the discretization of three- 
dimensional problems. Thanks to this simple shaped 
element, one can performed analytical derivations of the 
Cauchy principal value of the singular integrals. Such 
an approach has been previously applied for 2D elas- 
todynamic problems in (Dangla, 1988, 1990). To the 
author's best knowledge, this method has not been used 
in any 3D analysis. This paper addresses this issue by 
summarising the theoretical background of the method. 
Afterwards, the efficiency and accuracy of the regular- 
ization method is analysed in 3D elastodynamics. 



Numerical Modeling in Elastodynamics 

To analyze problems in 3D elastodynamics, various 
numerical methods are available: 

• the finite element method which is efficient to deal 
with complex geometries and numerous heterogene- 
ities (Chammas et al, 2003), even for inelastic con- 
stitutive models (Bonilla, 2000). It has neverthe- 
less several drawbacks such as numerical disper- 
sion (and damping) (Ihlenburg and Babuska, 1995; 
Semblat and Brioist, 2000) and (consequently) nu- 
merical cost in 3D elastodynamics, 

• the finite difference method which is very accurate 
in elastodynamics but is mainly adapted to simple 
geometries and linear constitutive models (Frankel 
and Vidale, 1992; Moczo et al, 2002; Virieux, 
1986) 

• the boundary element method which allows a very 
good description of the radiation conditions but is 
preferably dedicated to weak heterogeneities and 
linear constitutive models (Banerjee et al, 1988; 
Beskos, 1997; Beskos et al, 1986; Bonnet, 1999; 
Dangla, 1988; Sanchez-Sesma and Luzon, 1995; 
Yokoi, 2003) 

• the spectral element method which has been in- 
creasingly considered to analyse 2D / 3D wave prop- 
agation in linear media (Faccioli et al, 1996; Ko- 
matitsch and Vilotte, 1998) 

• the Aki-Larner method which takes advantage of 
the frequency-wavenumber decomposition but is 
limited to simple geometries (Aki and Larner, 1970; 
Bouchon et al, 1989) 

• series expansions of wave functions which give a 
semi-analytical estimation of the scattered wave- 
field for simple geometries (Moeen-Vaziri and Tri- 
funac, 1985; Sanchez-Sesma, 1983) 

Each method has specific advantages and drawbacks. It 
is consequently often more interesting to combine two 
methods to take advantage of their peculiarities. One of 
the most common method in elastodynamics is to couple 
FEM and BEM allowing an accurate description of the 
near field (FEM model including complex geometries, 
heterogeneities and constitutive behaviours) and a reli- 
able estimation of the far-field (BEM model involving 
radiation conditions). 

Integral Equations 

This paper is limited to isotropic elastodynamics for 
time-harmonic problems of circular frequency co. For 



any given body force distribution J^(x) over il, the gov- 
erning equations which must be verified by any displace- 
ment and stress fields, itj(x) and <xy (x), take the follow- 
ing form: 



aijj + puP'Ut + Fi = (2) 

The fundamental solutions, in time-harmonic elastody- 
namics, are defined by a force of unit amplitude applied 
at a fixed point y and in a fixed coordinate direction fc: 
Fj(x) = S(x — y)5ik- For infinite body the fundamen- 
tal solution, denoted by u^(x) = [/^'(x, y; a;), is known 
as the Helmholtz fundamental solution and is given by 
(Eringen and Suhubi, 1975): 



C/f(x,y;w) = 
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(3) 



where r 2 = (x — y) 2 and where kp = uj^J p/(\ + 2/i) 
and ks = oj^J p/fi are the longitudinal and transver- 
sal wave numbers respectively. The stress tensor asso- 
ciated with £/j fe (x, y; u>), defined by (TjQ), is denoted by 

(x, y; while the stress vector applied to the sur- 
face boundary of is T*(x, y; uj) = (x, y; uj)nj. 

For sake of simplicity let us assume no body force 
from now on. Application of the Maxwell-Betti reci- 
procity theorem leads to the following displacement in- 
tegral representation at point y € 1Z 3 (Bonnet, 1999): 



«(y)«*(y) = / [*i(x)0?(x,y;w) 

Jdfl 

-Ui(x)Tf(x,y;u;)] dS x 



(4) 



where k = 1 (y € 0) or k = (y ^ fi). 

Let y denote a fixed point on the boundary surface 
dfl. For a given small e > 0, introduce a spherical 
shaped neighbourhood u e (y) of y, called an exclusion 
neighbourhood (Guiggiani et ah, 1992). The domain 
^e(y) = ^ — v s(y) obtained by removing v e (y) from Q is 
such that the point y is exterior to f2 e (y) . Its boundary 



The free term (y) appearing in ©, is defined by: 

Cf(y) = lim/ ^(x,y;u,)dS x (7) 



(1) It is found to be equal to l/2S ik when f2 is smooth at y. 



Discretization and Regularization Principle 

The boundary surface dfl and the boundary vari- 
ables (ui,ti) are discretized by using three-noded trian- 
gular elements. A finite set of equations is generated 
by enforcing equation © at the nodes of the surface 
mesh (collocation method). Thus the boundary surface 
consists of the set of N boundary surface elements E e : 
dfl = {E e .e — 1..JV}. The integral appearing in © 
then assumes the form of a sum of TV element integrals: 



/ <■>-£/ 

Jdn „-tJe 



(■) 



(8) 



The numerical evaluation of non singular element inte- 
grals that appear in ([8]) is usually based, like in finite el- 
ement methods, on Gaussian quadrature formulas. The 
approximate value of an element integral can be given 
formally by: 



/(x) dS x 



/(x) = (9) 

=i 



where x, and Wi are the coordinates and weights of the 

Nu 

Gauss points. The notation J stands for the numer- 
ical approximation of integrals. This special notation 
has been adopted to emphasize that in case of singular 

Nu 

integral -f ^ J , 

Since some of element integrals are singular, a straigh- 
forward evaluation of ([5]) based on Gaussian quadrature 
formulas will inevitably lead to some significant error. 
To correct this error, a new term R k (y) must be intro- 
duced in the numerical evaluation of ©: 



Nu 



is dn e = (dQ-e e )+s e , where e e = dflDv s , s e = Qndv e . Cf(y) Ul {y) = R k (y)+ f [ti(x)U^ (x,y; u) 

Thp rl^QQiral fnrm nf thp int.pirrnl Pfinat.inn rnnsiet.s in J dfl 



(10) 



The classical form of the integral equation consists in 
taking the limit e — > in the representation formula Q 
taken for the domain The limiting expression thus 
obtained is known as the Somigliana identity: 

^(yR(y)=/ [U(x)U?{x,y;w) (5) 
Jon 

-u i (x)T i fc (x,y;w)] dS x 

The notation -f stands for the Cauchy principal value 
of a singular integral, i.e. the limit: 



-it*(x)lf (x,y;a;)] dS x 

The regularization method proposed in this paper 
consists in deriving analytically the correction term by 
taking advantage of the simple shape of the triangle el- 
ements. To do so, let us introduce the Kelvin's funda- 
mental solution: 



167T/l(l — v) 



(r.ir,* + (3 - Av)5 lk ) (11) 



/ (-) = lim/ 

Jdn E_>u J(dn-e.) 



(6) 



and note (x, y) the stress tensor associated with so- 
lution J7 t - (x, y). It is noticed that the Helmholtz and 



Kelvin solutions have identical singularities: 

(x->y) f UH^y;co)-UH*,y) = 0(l) 
(Vw>0) \ E£(x,y;w)-E&(x,y)=0(l) 



(12) 



Thanks to this property the correction term R k (y) 
only needs to involve the Kelvin fundamental solutions. 
For a given point y £ dfl, introduce the index subset 
-^(y) = {e € [1;-W]>y € E>e} such that the integral over 
E e is singular for e € I(y) and non singular for e ^ 
/(y). Introduce df2 y = {E e ,e e -f(y)} the set of the 
neighborhood elements of y. Thus the correction term 
can be formulated in the following form: 



R k (y) 



(13) 
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It can be noticed that formulation lfT3|) is indepen- 
dent of the interpolation order since tj(y) and itj(y) only 
need to be evaluated at point y. Taking advantage of 
the simple shape of the three-noded triangle elements, 
we can derive analytical expressions of the singular in- 
tegrals appearing in (fT3|) , 



In particular, it can be shown that they are the sum 
of elementary contributions involving elements of I(y): 



iHy) = / T*(x,y)dS a = Yl iHr,E e ) (14) 

an y eG/(y) 

JKy)=j- u?(*,y)ds x = J2 Ji(y;E e ) (is) 

9n y eeJ(y) 

The analytical derivations of I k (y;E e ) and J k (y;E e ) 
are proposed in the appendix (equations (f30|l to (|35)l ). 

In a similar manner, the free term involves the Kelvin 
fundamental solutions and can be assumed in the form 
of a sum of free term elements involving elements of 

i(y)- 

Cf(y) = lim f Tf{*,y)dS x = £ C?(y;£? e ) (16) 

where the exact derivation of C k (y;E e ) is given in the 
appendix (equations ([22)1 to ([29]) ). 

Since the method of derivation of the correction term 
i? fe (x, y) is now established, the formulation lfT0|) can be 
considered as the regularized form of the initial integral 
equation |(5]). 

Numerical Implementation 



Both boundary and unknowns are discretized using 
three-noded flat triangles and interpolation techniques 
initially developed for the Finite Element Method. The 
discretization of the geometry and the unknowns is thus 
written, respectively, as follows (Bonnet, 1999): 



x(£)=X>*(S) x * o(x) 



(17) 



fc=i 



fe=i 



with x fe : the node coordinates, Nk'- the linear interpola- 
tion functions and a k : the nodal values of the displace- 
ment or traction unknowns. 

Thus, the set of scalar equations resulting from the dis- 
cretization of equations lfT0|) . enforced at the nodes of 
the mesh, has the following matrix structure: 



[A]{u} + [B]{t} = 



(18) 



where [A] and [B] are fully populated non symmetric 
matrices. {u},{t} are the "vectors" containing, respec- 
tively, the nodal values of Ui (y) and U (y) . The incorpo- 
ration of the boundary conditions consists in substitut- 
ing the prescribed nodal values of (ui,U) into {u}, {t} 
in Eq. lfl8|) . The columns of this matrix equation are 
reordered so as to have a matrix equation of the form: 



[K]{v} = {f} 



(19) 



where the vector {v} consists of the unknown compo- 
nents of {u}, {t}. The matrix [K] contains the columns 
of [A], [B] associated with those unknown components 
while the right-hand side {f } results from the multiplica- 
tion of the known components of {u}, {t} by the corre- 
sponding columns of the matrices [A], [B]. As shown in 
the following for unbounded media, the right hand side 
{f} can also involve a contribution due to an incident 
wavefield. The method has been implemented into the 
computer code CESAR-LCPC (Humbert et al, 2005) of 
the Laboratoire Central des Ponts et Chaussees (French 
Public Works Research Laboratory, Paris, France). 



Figure 1: Cavity under harmonic internal pressure: 
model description. 

Validation in Frequency-Domain 
Elastodynamics 

Example 1: Spherical cavity under harmonic internal 
pressure. 

Description of the problem and analytical solution. 
The first example (figure Q]) concerns a spherical cav- 
ity of radius R in a full elastic isotropic space under- 
going an internal harmonic pressure. The cavity mesh 
includes 320 triangular boundary elements (that is 162 
nodes) and a special generation process is considered 



to have a regular triangular mesh of the sphere start- 
ing from an icosahedron (Edouard et al., 1996) (also see 
next sections). Using the regularization method pro- 
posed herein, we have computed the displacement field 
around the cavity at various (normalized) frequencies. 



and corrected numerical solutions. The efficiency of the 
regularization method then appears very good since the 
direct computation of the singular integrals leads to very 
bad results. 



The validation of the numerical results is made con- 
sidering the analytical solution in terms of radial dis- 
placement u(r,uj) given by Eringen and Suhubi (1975) 
as follows: 



u(r, u) 



P(uj)R 3 (ikp - l/r).exp (ikp(r - R)) 
^r{l-ik P R-k 2 s R 2 /4) 



(20) 



where hp and ks are the longitudinal and transverse 
wavenumbers. 

This equation can be rewritten using normalized dis- 
tance x = r /R, normalized frequency r/p = kpR/ir (that 
is r}p = 2R/Ap, Ap being the longitudinal wavelength) 

and considering v(x, T)p) ~ 



P{uj)R 



. It leads to: 



(ittXVp ~ l)-exp{inT]p{x - 1)) 
(1 - im] P - Tr 2 ( 2 r/ P /4).x 2 



(21) 



with C = kp/k s = y/J^QM)J{2^2v) = l/y/3 (that is 
v=0,25). 



Figure 2: Normalized radial displacement i>(x,f?p) (real 
part) vs normalized distance x : comparison between nu- 
merical and analytical results for normalized frequencies 
tip =0,01, 0.50, 1.00 and 2.00. 

Comparisons between numerical and analytical re- 
sults. In figure [21 the real part of the normalized radial 
displacement v(x, Tjp) defined by equation lj2"Tj) is dis- 
played vs normalized distance x f° r both analytical and 
numerical solutions at normalized frequencies ?yp=0.01, 
0.50, 1.00 and 2.00. For the nearly static case (r? P =0.01) 
as well as the fully dynamic cases, the agreement be- 
tween the numerical results and the analytical ones is 
very good at all normalized distances. From this first 
simple example, the reliability and accuracy of the pro- 
posed method then appear very good. 

Efficiency of the regularization method. We will then 
investigate the efficiency of the regularization method 
itself by evaluating the correction term for the same me- 
chanical problem (figure [T]). A non regularized solution 
is computed by dropping the correction term in lfT0|). In 
figure O this non regularized solution is compared with 
both the regularized one and the analytical solution at 
normalized frequencies ?yp=0.50 (left) and 2.00 (right). 
These comparisons show that the numerical results with- 
out the analytical correction are far from both analytical 



Figure 3: Normalized radial displacement v(x,Vp) (real 
part) vs normalized distance x : numerical results with 
and without the analytical correction of the singular inte- 
grals for normalized frequencies r\p=0.50 (left) and 2.00 
(right). 

Example 2: Diffraction of a plane wave by a spherical 
cavity. 

Description of the problem and analytical solution. 
The second example deals with the diffraction of a plane 
P-wave (ui nc = Uq exp[i(kpx — cut)] with Uq=1), prop- 
agating along x axis, by a spherical cavity. The numeri- 
cal results are firstly computed in frequency domain and 
compared with analytical results. They are afterwards 
converted into time domain to characterize the scattered 
wavefield. 

As shown in figure IH we have computed the wave field 
around the cavity for various directions. The boundary 
element mesh of the cavity (2562 nodes) is generated 
the same way as in the previous case (Edouard et al., 
1996). This mesh has been refined since the wave field 
has much stronger variations compared to the previous 
example. In figured! the radial displacement u{r, 0, <j), oj) 
is displayed vs distance x = r /R f° r both analytical and 
numerical solutions at two different normalized frequen- 
cies rjp (rjp = 2R/Kp). Different azimuthes are also 
considered. The analytical solution in terms of radial 
displacement u r is given by Pao and Mow (1973) as well 
as Eringen and Suhubi (19750. 

Comparisons between numerical and analytical re- 
sults in frequency domain. The results are computed 
for various azimuthes (0i = (i — 1) x 45°, 1 < i < 5) 
and figure [4] displays the real part of the radial dis- 
placement vs normalized distance X — r /R(^ — X — 



2 There are two mistakes in the original book of Eringen and 
Suhubi (1975) which have to be corrected as follows. The Origi- 
ns} 

nal expression of T X1 in (Eringen and Suhubi, 1975) (page 914, 
Eq. (9.12.11)) is: 

T x ( f (or) = {n 2 - n - \P 2 r 2 ) h£ ] (ar) + 2arh^ ] (or) 

and should be replaced by the following expression: 

Tj ( f(ar) = (n 2 - n - §/3V) (ar) + iarh^ (ar) 

The original expression of lC n in (Eringen and Suhubi, 1975) (page 
914, Eq. (9.12.13)) is: 

lC n = (l/A n )<^oi n (2n + l) \^$(aa)T$(aa)-T$(aa)T$(aa)\ 
and should be replaced by: 

lC n = (l/A n )</> i n (2n + l) \T^ ) ( a a)Tlf(aa)-T^ ) (aa)T[f (ao)] 
where a is the cavity radius, denoted R in this paper. 



Figure 4: Diffraction of a plane wave by a spheri- 
cal cavity: comparison with analytical results for vari- 
ous azimuthes at normalized frequencies r\p=t.00 and 
r]p =2.00. 

3) at two different normalized frequencies 77p=1.00 and 
?7p=2.00. The analytical results are plotted with lines 
(dotted for ?7p = 1.00 and solid for r?p=2.00) and the nu- 
merical results with symbols (circles for ?7p=1.00 and 
bullets for ?yp=2.00). The agreement between the nu- 
merical and analytical results is very good for all az- 
imuthes at ?7p=1.00. For 6*3 = 90°, some slight differ- 
ences can be noticed at ?7p=2.00 near the cavity wall. 
This is probably due to the fact that there is a grazing 
incidence at this point. 

Scattered wavefield in time domain. The numerical 
solutions are then estimated for various frequencies to 
compute the time domain scattered wavefield around the 
spherical cavity. As shown in figure El a Ricker signal, 
with normalized frequency ?yp=0.50, is considered for 
the excitation in time domain and the results are dis- 
played for three different azimuthes: = 0° correspond- 
ing to the direction of propagation (x axis), = 45° 
and = 90° (y axis) that is perpendicular to the direc- 
tion of propagation. For each azimuth, the time domain 
results are displayed for the upstream part of the prop- 
agation (x < —1) and the downstream part (\ > +!)• 
For the incident and transmitted wavefields, the vari- 
ous azimuthes do not always coincide with the direction 
of propagation whereas they correspond to directions of 
propagation of the scattered wavefield (see following ex- 
planations). As shown in figure the characterization 
of the scattered wavefield can then be easily performed 
as follows: 

• = 0° (top left): for this azimuth, only the X 
component of the displacement is displayed since 
the (computed) Y component is found negligible. 
The backward and forward components of the scat- 
tered wavefield clearly appear in the figure. For 
the upstream part, the scattered wavefield com- 
prises a P-wave as well as a S-wave component of 
respective velocities very close to the theoretical 
values (a few For the downstream part, the 
transmitted wavefield is easily identified and the 
scattered S-wave component has a velocity close to 
the previous value. Nevertheless, the amplitudes 
of the scattered wavefield components are not so 
large to identify them from figure [5j 

• 9 = 45° (center): for this azimuth, the apparent 
velocity of the incident and transmitted P-waves 

3 despite the fact we have considered a less refined mesh than 
for frequency domain computations. 



is lower because it does not coincide with the di- 
rection of propagation. Whereas for the scattered 
wavefield, radial directions correspond to the di- 
rection of propagation and the time domain nu- 
merical results show a large amplitude for both 
X and Y components. For the Y component of 
the scattered wavefield, both P and S-wave com- 
ponents can be identified in figured! The veloc- 
ity values estimated from the numerical results are 
found very close to theoretical ones. The velocity 
discrepancy between the downstream S component 
of the scattered wavefield and the transmitted P- 
wave is only due to the change of the apparent 
velocity of the latter which is azimuth dependent. 

• = 90° (bottom): for this azimuth, the apparent 
velocity of the incident P-wave is zero because it 
is perpendicular to the direction of propagation. 
The X and Y components of the displacement are 
displayed on one side of the cavity only since they 
are symmetrical on the opposite side. The X com- 
ponent clearly shows the S-wave part of the scat- 
tered wavefield. The estimation of its velocity is as 
good as in previous cases. For this azimuth, the Y 
component shows that the interaction between the 
plane wave and the cavity is particularly complex 
since we have a grazing incidence on the cavity 
wall. 



Figure 5: Diffraction of a plane P-wave by a spherical 
cavity: numerical results in time domain (Ricker signal) 
for various azimuthes at normalized frequency r]p=0.50. 

Example 3: Amplification of a plane seismic wave by 
a semi-spherical alluvial basin. 

Description of the problem and reference solution. 
The third example investigates the amplification of a 
plane seismic wave in an alluvial basin. In seismology 
and earthquake engineering, this phenomenon is known 
as "site effects" and generally leads to a strong ampli- 
fication of the seismic motion in soft alluvial deposits 
(Bard and Bouchon, 1985; Bielak et al., 1999; Chavez- 
Garcia et al, 2000; Moeen-Vaziri and Trifunac, 1985; 
Semblat et al, 2000, 2003a, 2005). The example con- 
sidered herein corresponds to a semi spherical alluvial 
basin (that is a soft elastic inclusion) in an elastic half 
space. Numerous papers have investigated the 3D wave 
diffraction by a semi spherical canyon (Lee, 1978; Liao et 
al, 2004; Yokoi, 2003) or 3D seismic wave amplification 
by surface heterogeneities (Dravinski, 2003; Komatitsch 
and Vilotte, 1998; Moczo et al, 2002; Sanchez-Sesma, 
1983; Sanchez-Sesma and Luzon, 1995). 

Several results have been published for the case of 
a semi spherical alluvial basin (Dravinski, 2003; Lee, 



Figure 6: Amplification of a plane vertical P-wave by a 
semi- spherical basin: model description. 

1984; Sanchez-Sesma, 1983). The 3D BEM model con- 
sidered herein for purpose of validation is depicted in 
figure [6j the mesh includes the semi-spherical basin of 
radius R (same type of triangular meshing as in the pre- 
vious section (Edouard et al, 1996)) and part of the 
free-surface (for r < 5R) . The contribution of the free 
surface r > 5i? in the BIE is neglected. Therefore, the 
BIE are enforced at the nodes of the mesh except those 
located at its boundary. The model is excited by a ver- 
tical plane P-wave. For the comparison, we will consider 
the results of Sanchez-Sesma (1983) derived thanks to a 
series expansion method. We will then investigate the 
amplification of the motion at the surface of the alluvial 
basin (i.e. soft inclusion). 

For the semi-spherical basin and the half-space, the 
mechanical parameters are chosen identical to Sanchez- 
Sesma's values as follows: 

• shear moduli: /j,r/ he = 0.3 

• mass densities: pr/pe = 0.6 

• Poisson's ratios: vr — 0.30 and ve — 0.25 

where subscript R refers to the alluvial basin and sub- 
script E to the half-space. 

Similarly to the previous section, we consider for the 
computations the same normalized frequency as Sanchez- 
Sesma corresponding to the diameter-to- wavelength ra- 
tio rjp = 2R/Ap where Ap is the P wavelength in the 
alluvial basin. 

Comparison between numerical and reference results. 
In figure [7j the amplification of the seismic motion is 
computed at the free-surface (vertical displacement) and 
displayed vs normalized distance (0 < % < 3). It is 
compared with Sanchez-Sesma's results (1983) at nor- 
malized frequency rjp = 0.50. The amplification of the 
vertical motion at the center of the semi-spherical basin 
is very well estimated by our numerical approach: 2.81 
for our numerical approach (i.e. 5.63 in amplitude to be 
compared to 2 for the half-space) and 2.82 for Sanchez- 
Sesma's results (i.e. 5.64 in amplitude). The com- 
puted displacement/distance curve from our numerical 
approach is very close to Sanchez-Sesma's semi-analytical 
results (figure [7]). This amplification value is larger than 
for the constant depth layer case (ID) since, for the semi- 
spherical basin, focusing effects are very strong (Sanchez- 
Sesma, 1983; Semblat et al, 2000, 2005). 
It should be noticed that the normalized frequency rjp = 
0.50 corresponds to the fundamental frequency of the ID 
case (the wavelength being Ap = 4R with R the depth 



of the basin). Nevertheless, at this frequency, the vari- 
ation of the amplification factor vs frequency is strong: 
for the 3D semi-spherical basin, this frequency is rather 
far from the maximum amplification peak. If we com- 
pute the amplification factor at the centre of the semi- 
spherical basin for various frequencies, the largest site 
effects are found at normalized frequency rjp = 0.57. At 
this frequency and for the mechanical properties chosen 
herein, the corresponding amplification factor is about 
4.76 (i.e. 9.52 in amplitude), that is 70% larger than for 
r)p = 0.50. For sake of comparisons, around normalized 
frequency rj P = 0.57 the amplitude variation with fre- 
quency is smaller (resp. basin properties). 



Figure 7: Computed vertical motion showing the amplifi- 
cation at the basin surface and comparison with Sanchez- 
Sesma's result for normalized frequency rjp — 0.50. 



Conclusion 

In this paper, a simple and efficient method to reg- 
ularize singular integrals in 3D boundary integral equa- 
tions has been presented. The regularization method is 
based on the derivation of analytical terms which correct 
the error due to the straightforward estimation of singu- 
lar integrals through classical Gaussian quadrature for- 
mulas. The analytical derivation of the correction term 
has assumed a three-noded triangle discretization of the 
boundary surfaces. However, the method described in 
the appendix can be easily generalized to any flat ele- 
ment such as quadrangle. This method has been im- 
plemented in a BEM code and applied to 3D frequency 
domain elastodynamics. 

Some comparisons have been made with (semi-)analy- 
tical results for simple problems: 

• cavity under harmonic pressure: the agreement 
between our numerical results and the analytical 
solution is very good for various frequencies even 
with a small number of nodes/elements. The ef- 
ficiency of the regularization method proposed in 
this paper is also discussed for this example. 

• diffraction of a plane-wave by a spherical cavity: 
the agreement between our numerical results and 
the analytical solution is very good for various az- 
imuthes and frequencies. In time domain, the nu- 
merical results are also found to be satisfactory. 

• wave amplification in a semispherical alluvial basin 
(soft inclusion): the comparison of our numeri- 
cal results with Sanchez-Sesma's semi-analytical 



results (Sanchez-Sesma, 1983) is also satisfactory. 
Further comparisons are planned with other cur- 
rent numerical approaches and more complex ge- 
ometries. 

Considering these good results, future work will then 
concern more realistic cases in the field of seismology. 
For sake of numerical efficiency, the regularization me- 
thod could also be implemented in a Symmetric Galerkin 
boundary element formulation (Bonnet et al, 1998) or in 
the framework of a Fast Multipole Method (Greengard 
et al, 1998). Our main goal is to have a detailed descrip- 
tion of the 3D geological structure of a given area to per- 
form reliable computations of seismic wave propagation 
and amplification (Bard and Bouchon, 1985; Bouchon 
et al, 1989; Chavez-Garcla et al., 2000; Frankel and 
Vidale, 1992; Moczo et al, 2002; Semblat et al, 2000, 
2003a,b, 2005). 



Appendix 
Calculation of C^(y;E e ) 

Let us calculate the free term Cf (y) defined by: 

C?(y) = lim f T?fry)dS a (22) 



In Eq. ([22]) . s £ is a spherical surface of radius e. Let 
x be a point on s e . The unit outward normal to s E is 
given by n = (y — x)/e. Thus the stress vector of the 
Kelvin fundamental solution applied to s £ has the form: 



?f(x,y) 



l 



l 



((1 - 2v)8 ik + 3n;n fe ) (23) 



8tt(1 - v) e 1 

Substituting this expression for T t k in ([22]) yields: 
l 



Cf(y) 



8tt(1-i/)' 



(24) 



lim 



(l — 2i/)\s e \S ik +3 / mn k dS. 



where \s £ \ is the surface area of s e . Here |s s | = s 2 ip, 
where ifi is the solid angle. A small amount of calcula- 
tions allows to derive the following expression: 



ni n k dS x = ^-5 ik 



r 2 E 

ee/(y) 



sin 



(25) 



where 9 e is the angle formed by the edges of E e at y, 
b\ is the unit vector of the bissecting line and is the 
unit outward normal to E e . The symmetry with respect 
to subscripts i and k in Eq. lf25|) shows the following 
identity: 



E 

ee/(y) 



e£/(y) 



sin — 6' 



(26) 



Combining ([24]) and ([25]) yields: 



Cf(y) = ls ik 



(27) 



Finally the solid angle is assumed to be the sum of ele- 
ment solid angles: 



(28) 



eeJ(y) 



Pratically, as shown in figure[Hl - e can be defined by the 
solid angle of the trihedron of apex y formed by the two 
edges of element E e and the semi-axis in the direction of 
— n(y), where n(y) is an arbitrary outward unit vector 
at y. In this case ip e = (if>\ + tp% + ip\ — 7r) where the ip\ 
are the three angles formed by the plane of the trihedron 
(figure [9]) . The calculation of the solid angles ip e relies 



on the knowledge of an outward unit vector at each node 
of the mesh. Practically for each node of coordinate y, 
n(y) can be calculated as the mean of unit normals to 
each element of J(y ) . It should be noticed that the aver- 
aging of the normal is only conventional. It results from 
an arbitrary choice in order to perform the calculation 
of solid angles ip e . The accuracy of the method does not 
depend on this averaging procedure since the value of 
the solid angle ip (equation l|28p ) is eventually recovered 
whatever the choice of n(y). Therefore, the free term 
Cj fe (y) is really the sum of free term elements C k (y; E e ) 
of the form: 



1 



8tt(1 - v) 



(29) 

(btn% + b%nt) 



Figure 8: Conical surface with apex at y 



then the sum of element integrals defined by: 

lHr,E e )= (33) 

— — -/ (nle k {a)-n%ei{a))\nL{a)da 

8tt(1 - v) J 



Figure 10: Distance L(a) on E e 
Calculation of J k (y; E e ) 

The integral J,- (y) can be written in the form: 

J*(y) = limV f U?(x,y)dS x (34) 

J(E, -e°) 



eel(y) 



Substituting expression (HU) for U?(x,y) in ([34]) gives 
the expression of J* (y; E e ): 



1 



(35) 



167T^(1 — v) Jq 



(eie k + (3 - Av)5 ik )L(a)da 



Figure 9: Description of angles tpi, ip% and if3 defining 
the element solid angle ip e . 

Calculation of ifXy; E e ) 

The integral /*(y) can be written in the form: 

I*(y) = limV / T t k (x,y)dS x (30) 



eel(y) 



where e e e = E e n v £ (with of course Yleeity) e t = e e)- 
Given e > 0, let us calculate the element integral ap- 
pearing in (|30| . Let x be a current point on E e and 
note n\ the unit outward normal to E e . The stress vec- 
tor of the Kelvin fundamental solution applied to E e is 
given by: 

Ja ,yj 8tt(1-i/) r2 [6L) 

where e; = {yi — Xi)/r. A trivial integration of the above 
expression shows that: 



L 



n i ek ~ n k e i ^ 



(n^efc(a) — n k ei(a)) In L(a)da (32) 



-2sinl - 1 {ntbl-n%bX)\ne 

where L(a) is the length of the segment defined in the 
figure |(10]). Thanks to identity ([26]), integral I 4 fe (y) is 
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